rm(list=ls())

pacman::p_load(lubridate, tidyverse, haven, gtrendsR, rdrobust, rddtools,
               broom, ggpubr, gridExtra,readxl)

setwd('put_your_wd_here')

# read in data
df <- read_excel("datasets/flag_sales.xlsx")
df <- janitor::clean_names(df)
df <- df %>%
  select(date, description, quantity) 

# unique flags
flags <- unique(df$description)

# group and count
all <- df %>%
  group_by(date) %>%
  count()

# add zeros for missing data
dates <- tibble(
  date = seq(min(all$date), max(all$date), by='day')
)
all <- left_join(dates, all)
all$n[is.na(all$n)] <- 0
all$date <- ymd(all$date)

# running variable for each shooting
all$day_running <- as.numeric(all$date - ymd('2021-03-22')) # Atlanta + Boulder
all$day_running_sj <- as.numeric(all$date - ymd('2021-05-26')) # San Jose
all$day_running_oh <- as.numeric(all$date - ymd('2021-11-30')) # Oxford High

estout <- function(each,x){
  dat <- each
  sd_est <- sd(as.numeric(dat$n),na.rm=T)
  house_rdd <- rdd_data(y=as.numeric(dat$n), x=pull(dat[,x]), cutpoint=0)
  bw_ik <- rdd_bw_ik(house_rdd)
  reg_nonpara <- rdd_reg_np(rdd_object=house_rdd, bw=bw_ik)
  rddtoolsest <- data.frame(Est = 'Rdd Tools Coef',
                            Coef=reg_nonpara$coefficients/sd_est,
                            lower=(reg_nonpara$coefficients - 1.96 * reg_nonpara$coefMat[2])/sd_est,
                            upper=(reg_nonpara$coefficients + 1.96 * reg_nonpara$coefMat[2])/sd_est)
  return(rddtoolsest)
}

# each flag separately
each <- df %>%
  group_by(date, description) %>%
  count()
each <- left_join(dates, each)
each$n[is.na(each$n)] <- 0
each$date <- ymd(each$date)
each$day_running <- as.numeric(each$date - ymd('2021-03-22')) # atlanta 
each$day_running_sj <- as.numeric(each$date - ymd('2021-05-26')) # san jose
each$day_running_oh <- as.numeric(each$date - ymd('2021-11-30')) # Oxford HIgh

plot1 <- bind_rows(
  estout(each = all %>% filter(day_running >= -60 & day_running <= 60),'day_running'),
  estout(each = each %>% filter(day_running >= -60 & day_running <= 60 & description == flags[1]),'day_running'),
  estout(each = each %>% filter(day_running >= -60 & day_running <= 60 & description == flags[2]),'day_running'),
  estout(each = each %>% filter(day_running >= -60 & day_running <= 60 & description == flags[3]),'day_running')
) %>%
  mutate(outcome = c('All Flags','\"Come and Take It\"', '\"Liberty or Death\"', '\"2nd Amendment\"'),
         outcome = factor(outcome, levels=rev(c(
           c('All Flags','\"Come and Take It\"', '\"Liberty or Death\"', '\"2nd Amendment\"')
         )))) %>%
  ggplot(aes(x=outcome, y=Coef, ymin=lower, ymax=upper)) +
  geom_errorbar(width=0.01) +
  geom_point(shape=21, size=3, fill='white') +
  geom_hline(yintercept=0, color='red', linetype=2) +
  coord_flip() + 
  theme_bw() +
  labs(x='', y='Effect on Flag Sales\n(Standard Deviations)') +
  theme(axis.title.x = element_text(margin=margin(t = 10))) +
  scale_y_continuous(limits=c(-2, 2), breaks=seq(-2,2,by=0.5))
plot1
ggsave(plot1, 
       width=4, height=2.5, 
       filename = 'figures/flagsalesRDIT.png')

# other shootings
out1 <- bind_rows(
  estout( each = all %>% filter(day_running_sj >= -60 & day_running_sj <= 60),'day_running_sj'),
  estout( each = each %>% filter(day_running_sj >= -60 & day_running_sj <= 60 & description == flags[1]),'day_running_sj'),
  estout( each = each %>% filter(day_running_sj >= -60 & day_running_sj <= 60 & description == flags[2]),'day_running_sj'),
  estout( each = each %>% filter(day_running_sj >= -60 & day_running_sj <= 60 & description == flags[3]),'day_running_sj')
) %>%
  mutate(outcome = c('All Flags','\"Come and Take It\"', '\"Liberty or Death\"', '\"2nd Amendment\"'),
         outcome = factor(outcome, levels=rev(c(
           c('All Flags','\"Come and Take It\"', '\"Liberty or Death\"', '\"2nd Amendment\"')
         ))),
         Shooting = 'San Jose VTA') 
out2 <- bind_rows(
  estout( each = all %>% filter(day_running_oh >= -60 & day_running_oh <= 60),'day_running_oh'),
  estout( each = each %>% filter(day_running_oh >= -60 & day_running_oh <= 60 & description == flags[1]),'day_running_oh'),
  estout( each = each %>% filter(day_running_oh >= -60 & day_running_oh <= 60 & description == flags[2]),'day_running_oh'),
  estout( each = each %>% filter(day_running_oh >= -60 & day_running_oh <= 60 & description == flags[3]),'day_running_oh')
) %>%
  mutate(outcome = c('All Flags','\"Come and Take It\"', '\"Liberty or Death\"', '\"2nd Amendment\"'),
         outcome = factor(outcome, levels=rev(c(
           c('All Flags','\"Come and Take It\"', '\"Liberty or Death\"', '\"2nd Amendment\"')
         ))),
         Shooting = 'Oxford High') 

out3 <- bind_rows(
  estout(each = all %>% filter(day_running >= -60 & day_running <= 60),'day_running'),
  estout(each = each %>% filter(day_running >= -60 & day_running <= 60 & description == flags[1]),'day_running'),
  estout(each = each %>% filter(day_running >= -60 & day_running <= 60 & description == flags[2]),'day_running'),
  estout(each = each %>% filter(day_running >= -60 & day_running <= 60 & description == flags[3]),'day_running')
) %>%
  mutate(outcome = c('All Flags','\"Come and Take It\"', '\"Liberty or Death\"', '\"2nd Amendment\"'),
         outcome = factor(outcome, levels=rev(c(
           c('All Flags','\"Come and Take It\"', '\"Liberty or Death\"', '\"2nd Amendment\"')
         ))),
         Shooting = 'Atlanta + Boulder') 
dat <- bind_rows(out1, out2, out3) %>%
  filter(outcome != 'All Flags')
calcMeta <- function(data){
  m.gen <- meta::metagen(TE = Coef,
                         seTE = se,
                         studlab = Shooting,
                         data = data,
                         sm = "SMD",
                         fixed = FALSE,
                         random = TRUE)
  results <- tibble(Coef = m.gen$TE.random,
                    lower = m.gen$lower.random,
                    upper = m.gen$upper.random)
  return(results)
}
meta1 <- calcMeta(dat %>% mutate(se = (upper - Coef)/1.96)) %>%
  mutate(Shooting='Meta Estimate',
         outcome='Meta Analysis Estimate')
dat <- bind_rows(dat, meta1)
plot2 <- dat %>%
  mutate(statsig = ifelse((lower > 0 & upper > 0) | (lower < 0 & upper < 0),1,0)) %>%
  mutate(Shooting = factor(Shooting, levels=c('Meta Estimate','Oxford High', 'San Jose VTA','Atlanta + Boulder'
                                              ))) %>%
  mutate(outcome = factor(outcome, levels=c('Meta Analysis Estimate','\"Come and Take It\"', '\"Liberty or Death\"', '\"2nd Amendment\"'))) %>%
  ggplot(aes(x=Shooting, y=Coef, ymin=lower, ymax=upper, shape=outcome, color=factor(statsig))) +
  geom_hline(yintercept=0, color='red', linetype=2) +
  geom_errorbar(width=0.01, position=position_dodge(width=0.3)) +
  geom_point(size=2, fill='white', position=position_dodge(width=0.3)) +
  coord_flip() + 
  theme_bw() +
  labs(x='', y='Effect on Flag Sales\n(Standard Deviations)') +
  theme(axis.title.x = element_text(margin=margin(t = 10))) +
  scale_y_continuous(limits=c(-3, 3), breaks=seq(-3,3,by=1)) +
  scale_shape_manual(values=c(21, 23, 24, 16), name='Flag Slogan',
                     breaks=rev(c('Meta Analysis Estimate','\"Come and Take It\"', '\"Liberty or Death\"', '\"2nd Amendment\"'))) +
  scale_color_manual(values = c('grey50','black'), guide='none')   
  
plot2

ggsave(plot2, 
       width=6, height=3.5, 
       filename = 'figures/flagsalesRDIT_extra.png')

# meta analysis
datformeta <- bind_rows(out1, out2, out3, out4) %>%
  filter(outcome != 'All Flags') 
datformeta$se = (datformeta$upper - datformeta$Coef)/1.96

m.gen <- meta::metagen(TE = Coef,
                 seTE = se,
                 studlab = Shooting,
                 data = datformeta,
                 sm = "SMD",
                 fixed = FALSE,
                 random = TRUE)
results <- tibble(pe = m.gen$TE.random,
                  lower = m.gen$lower.random,
                  upper = m.gen$upper.random)

